%% load calibration South at initial steady state
str_region="South_";
str1="params_";
 str_3="calibration_Table2";
 str_save=append(path_calibration,str1,str_region,str_3,".mat");
filename=cell2mat(str_save);
load(filename)
read_pars;read_pars_s;
n_p=length(pars);
pars(n_p+1)=incid; pars(n_p+2)=k;
pars(n_p+2)=k_s_ben;%k_s_ben=k_s;
A=A_s_ben;

Z_s=z_0_s*reshape(exp(z_vec)'*ones(1,N_b),Ntot,1);

%% compute nwe steady state and dynamics
welf_disc_all=[];
ip2=0;inde_solved = ones(length(tau_vec),length(k_vec));
while ip2<length(k_vec)
    ip2
    ip2=ip2+1;
    size_subsidy=1-k_vec(ip2);
    incid = incid_vec(ip2);
    ip=0;
    while ip<length(tau_vec)
        ip=ip+1;
        ip
        global cub_spl_v cub_spl_x pp_star Amin

        %find the equilibrium profitability of firms that is consistent
        %withthe entry cost
        k_s = (1-size_subsidy*(1-tau_vec(ip)))*k_s_ben;
        tau_s = tau_vec(ip)*size_subsidy*k_s_ben;%part ex-post
        pars(7)=k_s;pars(10)=tau_s;pars(n_p+1)=incid;pars_s=pars;
        
        %% find new value added and debt at entry
        n_try=0;
        fv=ones(3,1);A_0=A_s_ben*0.99; if tau_vec(ip)>0.4,B_0=1.2*B0_s;else,B_0=0.25*B0_s;end
        if ip==1,vars_0=[B_0,A_0,B_0];end
        while max(abs(fv))>1/10000 && n_try<20
            n_try=n_try+1;
            [vars,fv]=fsolve('fun_EntryEq',vars_0,options,pars_s);
            if incid==1
                A_0=A_0*0.95;B_0=B0_s;
                vars_0=[B_0,A_0,B_0];
            else
                vars_0=[0.8*B_0,A_0,1.2*B_0];
            end
        end

        B0_s     = max(0,vars(1));A_s_new=vars(2);B00_s=max(0,vars(3));
     
        b_common=B0_s/(z_0*A_s_new);
        b0_low_s=b_common/zhig_s;
        b0_mid_s=b_common/zmid_s;
        b0_hig_s=b_common/zlow_s;
        %compute debt buyback when rebate is expost
        if tau_s>0 && b0_hig_s>ppval(pp_star,A_s_new)
            b0_hig_s=fun_buyback(b0_hig_s,A_s_new,tau_s/(z_0*A_s_new),cub_spl_v,cub_spl_x);
        end
        
        if max(abs(fv))>1/1000 || b0_hig_s>ppval(pp_bar,A_s_new)
            inde_solved(ip,ip2)=0;%this is an indicator to look at ex-post. it's equal to 1 if an equilibrium exists
        end

        if tau_s==0 && size_subsidy==0,b_up_ss=b_up; b_md_ss=b_md;b_dn_ss=b_dn;end
        b_up=b0_hig_s;b_md=b0_mid_s;b_dn=b0_low_s;
        %untreated firms
        b_common0 = B00_s/(z_0*A_s_new);
        b_dn2=b_common0/zhig_s;
        b_md2=b_common0/zmid_s;
        b_up2=b_common0/zlow_s;

        %% find invariant distribution (b,z) at new steady state in the South
        b_bar_s=ppval(pp_bar,A_s_new);if ip==1 && ip2==1,b_bar_ss=b_bar_s;end
        q_hig=q_hig_s;q_low=q_low_s;q_mid=q_mid_s;
        bhig=b_up;bmid=b_md;blow=b_dn;
        b_bar=b_bar_s;
        get_distribution_b;
        inde_med = find(cumsum(f_dist)/sum(f_dist)>0.5,1)-1;
        med_debt_coupon_s(ip,ip2) = (coupon)./ppval(pp_x,b_vec_d(inde_med));
        b_vec_mat_l_s = b_vec_mat_l+log(b_bar_s);
        b_vec_mat_s   = exp(b_vec_mat_l_s);
        b_mat_s       = ones(length(z_vec),1)*(b_vec_mat_s);
        bb_vec_mat_s  = reshape(b_mat_s,N_z*N_b,1);
        b_s           = reshape(ones(N_z,1)*(b_vec_mat_s),Ntot,1);
        inv_s         = zeros(Ntot,1);
        for ii=1:Ntot,inv_s(ii) = fnval(cub_spl_i,[b_s(ii);A_s_new]);end
        b_vec_mat=b_vec_mat_s;bb_vec_mat=bb_vec_mat_s;
        get_distribution_bZ;
        f_dist_s          = (f_dist);
        P_default_avg_s   = 1/sum(f_dist_s)+delta_s;
        debt_s            = (bb_vec_mat'*f_dist_s)/sum(f_dist_s);
        A_T_s   = A_T;
        h_00_s  = h_00;
        inde_b_bar_s=inde_b_bar;
        inde_bi_s=inde_bi;
        inde_bb_s=inde_bb;
        inde_nego_s= inde_nego;
        if ip==1 && ip2==1
            A_T_s_1   =A_T;
            h_00_s_1  =h_00;
            sun_s_1=sum(f_dist);
        end
        avg_debt_coupon_s(ip,ip2) = (((coupon)./ppval(pp_x,bb_vec_mat))'*f_dist_n)/sum(f_dist);

        %% new aggregate steady state
        Y_new = (wei_s*l_s_ss*A_s_ben^(1/(1-nu)) + wei_n*l_n_ss*A_n_ben^(1/(1-nu)))^((nu-1)/(nu-2));
        crit_Y=1;l_n=l_n_ss;l_s=l_s_ss;
        while crit_Y>1/100000
            w_s       = (nu-1)/nu*(A_s_new/Y_new)^(1/(1-nu));
            w_n       = (nu-1)/nu*(A_n_new/Y_new)^(1/(1-nu));
            if l_reallocation==1
                [l_n,fv]      =  fsolve(@(ell) w_n - w_s - h_n*max(0,ell).^eta  + h_s*max(0,((1-ell*wei_n)/wei_s)).^eta,1,options);
                if abs(fv)>1/10000,error('fv'),end
                l_s           = (1-l_n*wei_n)/wei_s;%l_s*wei_s+l_n*wei_n=1
            end
            if demand_ext==1
                Y_imp         =  (wei_s*l_s*A_s_new^(1/(1-nu)) + wei_n*l_n*A_n_new^(1/(1-nu)))^((nu-1)/(nu-2));
            else
                Y_imp         = Y_new;
            end
            crit_Y        = abs(Y_imp- Y_new) ;
            Y_new         = Y_new*(1+0.1*(Y_imp- Y_new)/Y_new);
        end
        %aggregate outcomes in each region
        m_s       = (A_s_new/Y_new)^(nu/(1-nu))/(Y_new*(sum(Z_s.*f_dist_s)))*sum(f_dist_s);
        m_n       = (A_n_new/Y_new)^(nu/(1-nu))/(Y_new*(sum(Z_n.*f_dist_n)))*sum(f_dist_n);
        f_s       = m_s*f_dist_s/sum(f_dist_s);
        f_n       = m_n*f_dist_n/sum(f_dist_n);
        m_tilde_s = m_s*P_default_avg_s;
        m_tilde_n = m_n*P_default_avg_n;
        inve_s    = A_s_new*sum(Z_s.*inv_s.*f_s);
        inve_n    = A_n_new*sum(Z_n.*inv_n.*f_n);
        y_s       = (A_s_new/Y_new)^(1/(1-nu));
        y_n       = (A_n_new/Y_new)^(1/(1-nu));
        i_s       =  k_s_ben*m_tilde_s;
        i_n       =  k_n_ben*m_tilde_n;

        Y_vec(ip,ip2)   = y_s*l_s*wei_s+y_n*l_n*wei_n;
        I_vec(ip,ip2)   = i_s*l_s*wei_s+i_n*l_n*wei_n;
        C_vec(ip,ip2)   = Y_vec(ip,ip2)-I_vec(ip,ip2);
        c_s_vec(ip,ip2) = y_s-i_s+i_s*size_subsidy*(1-sub_tax_s)-inve_s;
        c_n_vec(ip,ip2) = y_n-i_n-i_s*size_subsidy*(1-sub_tax_s)*l_s*wei_s/(l_n*wei_n)-inve_n;

        W_s_vec_ss(ip,ip2) = c_s_vec(ip,ip2)-chi_s*(sum(Z_s.*f_s))-h_s*l_s^(eta);
        W_n_vec_ss(ip,ip2) = c_n_vec(ip,ip2)-chi_n*(sum(Z_n.*f_n))-h_n*l_n^(eta);
        W_vec_ss(ip,ip2)   = W_s_vec_ss(ip,ip2)*l_s*wei_s+W_n_vec_ss(ip,ip2)*l_n*wei_n;

        bbar_vec_s(ip,ip2)= b_bar_s;
        bbar_vec_n(ip,ip2)= b_bar_n;
        R_vec_n(ip,ip2)= 1/nu*A_n_new-chi_n;
        R_vec_s(ip,ip2)= 1/nu*A_s_new-chi_s;
        A_vec_n(ip,ip2)= A_n_new;
        A_vec_s(ip,ip2)= A_s_new;
        B0_vec_s(ip,ip2)= B0_s;
        B0_vec_n(ip,ip2)= B0_n;
        l_s_vec(ip,ip2) = l_s;
        l_n_vec(ip,ip2) = l_n;
        dist_vec_n(:,ip,ip2) = f_n;
        supp_vec_n(:,ip,ip2) = bb_vec_mat_n/A_n_new;
        dist_vec_s(:,ip,ip2) = f_s;
        supp_vec_s(:,ip,ip2) = bb_vec_mat_s/A_n_new;
        y_s_vec(ip,ip2) = y_s;
        y_n_vec(ip,ip2) = y_n;
        inve_s_vec(ip,ip2) = (inve_s+i_s)/y_s;
        inve_n_vec(ip,ip2) = (inve_n+i_n)/y_n;
        Q_s_vec(ip,ip2) = ((w_s/(1-1/nu))^(1-nu))*(sum(Z_s.*f_s))*Y_new;
        Q_n_vec(ip,ip2) = ((w_n/(1-1/nu))^(1-nu))*(sum(Z_n.*f_n))*Y_new;
        M_s_vec(ip,ip2) = sum(f_s);
        M_n_vec(ip,ip2) = sum(f_n);
        m_s_vec(ip,ip2) = m_tilde_s;
        m_n_vec(ip,ip2) = m_tilde_n;
        va_s_vec(ip,ip2) = A_s_new;
        va_n_vec(ip,ip2) = A_n_new;
        z_s_vec(ip,ip2)= sum(Z_s.*f_s)/sum(f_s);
        z_n_vec(ip,ip2)= sum(Z_n.*f_n)/sum(f_n);
        w_s_vec(ip,ip2) = w_s;
        w_n_vec(ip,ip2) = w_n;
        P_exit_n_vec(ip,ip2) = P_default_avg_n;
        P_exit_s_vec(ip,ip2) = P_default_avg_s;
        debt_va_s_vec(ip,ip2)= debt_s;
        debt_va_n_vec(ip,ip2)= debt_n;
        debt0_va_s_vec(ip,ip2)= (zlow_s*b0_hig_s*q_hig_s+zhig_s*b0_low_s*q_low_s+zmid_s*b0_mid_s*q_mid_s)*A_s_new/R_vec_s(ip,ip2);
        debt0_va_n_vec(ip,ip2)= (zlow_n*b0_hig_n*q_hig_n+zhig_n*b0_low_n*q_low_n+zmid_n*b0_mid_n*q_mid_n)*A_n_new/R_vec_n(ip,ip2);
        debt0t_va_s_vec(ip,ip2)=(zlow_s*b_up2*q_hig_s+zhig_s*b_dn2*q_low_s+zmid_s*b_md2*q_mid_s)*A_s_new/R_vec_s(ip,ip2);

        if k_s==k_s_ben && tau_s==0
            f_dist_ss_n=f_n;f_dist_ss_s=f_s;
        end

        %% Obtain transition dynamics
        get_transition_dynamics;

    end

end
